# analysis_section_4_2_part1.R
# This script replicates the analysis in Section 4.2 (Part 1) of the paper
#   Figure 3, Figure 5, Table SI.5, Table SI.6, and Table SI.8
# Author: Yimeng Li
# Dependencies:
#   Data/VR_VoteCal_LA.Rdata

# Load libraries:
library(dplyr)
library(ggplot2)
library(patchwork)
library(sandwich)
library(tidyr)

# Load data:
load("./Replication Package/Data/VR_VoteCal_LA.Rdata")

# Convert data from wide to long format for analysis
VR_VoteCal_LA_long <- VR_VoteCal_LA %>%
  filter(RegisteredBefore2016Pri == 1, !is.na(Dist)) %>%
  select(Dist, UVBM, IsPermanentVBM, PoliticalParty, Turnout.2020, Turnout.2016, Vote_by_Mail.2020, Vote_by_Mail.2016) %>%
  pivot_longer(cols = c("Turnout.2020", "Turnout.2016", "Vote_by_Mail.2020", "Vote_by_Mail.2016"),
               names_to = c(".value", "Year"),
               names_sep = "[.]")

#*******************************************************************************
# Analysis - Define Functions
#*******************************************************************************

# This function runs the difference-in-differences analysis given
#   DV: Dependent variable (Vote_by_Mail, or Turnout)
#   Spec: Specification (Levels, Differences, or Diff-in-Diffs)
#   PVBM: Permanent VBM status (Non-Perm. VBM, Perm. VBM, or All Voters)
#   Party: Political party (Democratic, Republican, or NPA)
run_model <- function(DV, Spec, PVBM, Party) {
  
  IVs = case_when(
    Spec == "Levels" ~ '0 + interaction(UVBM, Year)',
    Spec == "Differences" ~ '0 + I(UVBM == 1) + I(UVBM == 0 & Year == "2020") + I(UVBM == 1 & Year == "2020")',
    Spec == "Diff-in-Diffs" ~ 'I(UVBM == 1) + I(Year == "2020") + I(UVBM == 1)*I(Year == "2020")'
  )
  
  formula = as.formula(paste(DV, "~", IVs))
  
  PVBM_set = case_when(
    PVBM == "Non-Perm. VBM" ~ c("Non-Perm. VBM"),
    PVBM == "Perm. VBM" ~ c("Perm. VBM"),
    PVBM == "All Voters" ~ c("Non-Perm. VBM", "Perm. VBM")
  )
  
  Party_set = case_when(
    Party == "Democratic" ~ c("Democratic"),
    Party == "Republican" ~ c("Republican"),
    Party == "NPA" ~ c("No Party Preference"),
    Party == "All Parties" ~ c("Democratic", "Republican", "No Party Preference", "Other")
  )
  
  data = VR_VoteCal_LA_long %>%
    filter(IsPermanentVBM %in% PVBM_set, PoliticalParty %in% Party_set)
  
  lm_model = lm(formula, data)
  
  return(lm_model)
}

# This function calculates the estimates and clustered standard errors given
#   lm_model: the difference-in-differences model output by run_model()
get_estimates <- function(lm_model){
  
  lm_estimates <- data.frame(
    term = attributes(lm_model$coefficients)$names,
    coef = as.numeric(lm_model$coefficients),
    se = sqrt(diag(vcovCL(lm_model, cluster = ~ Dist)))
  )
  
  return(lm_estimates)
}

# This function generates the difference-in-differences figures given
#   lm_estimates_NPAV: Estimates of Pre and Post levels for Non-Perm. VBM voters output by get_estimates()
#   lm_estimates_PAV: Estimates of Pre and Post levels for Perm. VBM voters output by get_estimates()
#   y_axis_limits: Y-axis limits
#   y_axis_breaks: Y-axis breaks
#   y_axis_text: Y-axis text
#   title_text: Title text
gen_figure <- function(lm_estimates_NPAV, lm_estimates_PAV, y_axis_limits, y_axis_breaks, y_axis_text, title_text) {
  
  lm_estimates_NPAV_PAV <- bind_rows(
    lm_estimates_NPAV %>%
      mutate(UVBM = c("No", "Yes", "No", "Yes"),
             IsPermanentVBM = c("No", "No", "No", "No"),
             Year = c("2016", "2016", "2020", "2020")),
    lm_estimates_PAV %>%
      mutate(UVBM = c("No", "Yes", "No", "Yes"),
             IsPermanentVBM = c("Yes", "Yes", "Yes", "Yes"),
             Year = c("2016", "2016", "2020", "2020"))
    ) %>%
    mutate(Group = case_when(IsPermanentVBM == "No" & UVBM == "No" ~ "Non-Perm. VBM Voters in Non-Univ. VBM Districts",
                             IsPermanentVBM == "No" & UVBM == "Yes" ~ "Non-Perm. VBM Voters in Univ. VBM Districts",
                             IsPermanentVBM == "Yes" & UVBM == "No" ~ "Perm. VBM Voters in Non-Univ. VBM Districts",
                             IsPermanentVBM == "Yes" & UVBM == "Yes" ~ "Perm. VBM Voters in Univ. VBM Districts"),
           Group = factor(Group, levels = c("Non-Perm. VBM Voters in Non-Univ. VBM Districts",
                                            "Perm. VBM Voters in Non-Univ. VBM Districts",
                                            "Non-Perm. VBM Voters in Univ. VBM Districts",
                                            "Perm. VBM Voters in Univ. VBM Districts")))
  
  Figure <- ggplot(lm_estimates_NPAV_PAV, aes(x = Year, y = coef, colour = Group, linetype = Group)) +
    geom_line(aes(group = Group), size = 1) +
    geom_point(size = 2, show.legend = FALSE) +
    scale_y_continuous(limits = y_axis_limits, breaks = y_axis_breaks, labels = scales::percent_format(accuracy = 1)) +
    scale_colour_manual(values = c("blue", "blue", "orange", "orange")) +
    scale_linetype_manual(values = c("dashed", "solid", "dashed", "solid")) +
    labs(x = "Year", y = y_axis_text, title = title_text, colour = "", linetype = "") +
    theme_bw() +
    theme(panel.border = element_rect(colour = "black", fill = NA, size = 1), panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
    theme(text = element_text(size = 14), axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
    theme(legend.position = "bottom", legend.direction = "vertical",
          legend.key.size =  unit(0.5, "in"), legend.text = element_text(size = 14)) +
    guides(colour = guide_legend(nrow = 2, byrow = TRUE))
  
  return(Figure)
}

#*******************************************************************************
# Difference-in-Differences Analysis without Covariates
# Results Reported in Figure 3, Table SI.5, and Table SI.6
#*******************************************************************************

# Estimate Effects on Percent Voting by Mail
# Non-Permanent VBM Voters
PercVBM_NPAV_Part1 <- run_model("Vote_by_Mail", "Levels", "Non-Perm. VBM", "All Parties")
PercVBM_NPAV_Part2 <- run_model("Vote_by_Mail", "Differences", "Non-Perm. VBM", "All Parties")
PercVBM_NPAV_Part3 <- run_model("Vote_by_Mail", "Diff-in-Diffs", "Non-Perm. VBM", "All Parties")

# Permanent VBM Voters
PercVBM_PAV_Part1 <- run_model("Vote_by_Mail", "Levels", "Perm. VBM", "All Parties")
PercVBM_PAV_Part2 <- run_model("Vote_by_Mail", "Differences", "Perm. VBM", "All Parties")
PercVBM_PAV_Part3 <- run_model("Vote_by_Mail", "Diff-in-Diffs", "Perm. VBM", "All Parties")

# Non-Perm. And Perm. VBM Voters
PercVBM_ALL_Part1 <- run_model("Vote_by_Mail", "Levels", "All Voters", "All Parties")
PercVBM_ALL_Part2 <- run_model("Vote_by_Mail", "Differences", "All Voters", "All Parties")
PercVBM_ALL_Part3 <- run_model("Vote_by_Mail", "Diff-in-Diffs", "All Voters", "All Parties")

# Obtain Estimates and Clustered Standard Errors
# Non-Permanent VBM Voters
PercVBM_NPAV_Part1_est <- get_estimates(PercVBM_NPAV_Part1)
PercVBM_NPAV_Part2_est <- get_estimates(PercVBM_NPAV_Part2)
PercVBM_NPAV_Part3_est <- get_estimates(PercVBM_NPAV_Part3)

# Permanent VBM Voters
PercVBM_PAV_Part1_est <- get_estimates(PercVBM_PAV_Part1)
PercVBM_PAV_Part2_est <- get_estimates(PercVBM_PAV_Part2)
PercVBM_PAV_Part3_est <- get_estimates(PercVBM_PAV_Part3)

# Non-Perm. And Perm. VBM Voters
PercVBM_ALL_Part1_est <- get_estimates(PercVBM_ALL_Part1)
PercVBM_ALL_Part2_est <- get_estimates(PercVBM_ALL_Part2)
PercVBM_ALL_Part3_est <- get_estimates(PercVBM_ALL_Part3)

# Number of Observations
N_PercVBM_NPAV <- nobs(PercVBM_NPAV_Part1)
N_PercVBM_PAV <- nobs(PercVBM_PAV_Part1)
N_PercVBM_ALL <- nobs(PercVBM_ALL_Part1)

# Estimate Effects on Turnout
# Non-Permanent VBM Voters
Turnout_NPAV_Part1 <- run_model("Turnout", "Levels", "Non-Perm. VBM", "All Parties")
Turnout_NPAV_Part2 <- run_model("Turnout", "Differences", "Non-Perm. VBM", "All Parties")
Turnout_NPAV_Part3 <- run_model("Turnout", "Diff-in-Diffs", "Non-Perm. VBM", "All Parties")

# Permanent VBM Voters
Turnout_PAV_Part1 <- run_model("Turnout", "Levels", "Perm. VBM", "All Parties")
Turnout_PAV_Part2 <- run_model("Turnout", "Differences", "Perm. VBM", "All Parties")
Turnout_PAV_Part3 <- run_model("Turnout", "Diff-in-Diffs", "Perm. VBM", "All Parties")

# Non-Perm. And Perm. VBM Voters
Turnout_ALL_Part1 <- run_model("Turnout", "Levels", "All Voters", "All Parties")
Turnout_ALL_Part2 <- run_model("Turnout", "Differences", "All Voters", "All Parties")
Turnout_ALL_Part3 <- run_model("Turnout", "Diff-in-Diffs", "All Voters", "All Parties")

# Obtain Estimates and Clustered Standard Errors
# Non-Permanent VBM Voters
Turnout_NPAV_Part1_est <- get_estimates(Turnout_NPAV_Part1)
Turnout_NPAV_Part2_est <- get_estimates(Turnout_NPAV_Part2)
Turnout_NPAV_Part3_est <- get_estimates(Turnout_NPAV_Part3)

# Permanent VBM Voters
Turnout_PAV_Part1_est <- get_estimates(Turnout_PAV_Part1)
Turnout_PAV_Part2_est <- get_estimates(Turnout_PAV_Part2)
Turnout_PAV_Part3_est <- get_estimates(Turnout_PAV_Part3)

# Non-Perm. And Perm. VBM Voters
Turnout_ALL_Part1_est <- get_estimates(Turnout_ALL_Part1)
Turnout_ALL_Part2_est <- get_estimates(Turnout_ALL_Part2)
Turnout_ALL_Part3_est <- get_estimates(Turnout_ALL_Part3)

# Number of Observations
N_Turnout_NPAV <- nobs(Turnout_NPAV_Part1)
N_Turnout_PAV <- nobs(Turnout_PAV_Part1)
N_Turnout_ALL <- nobs(Turnout_ALL_Part1)

#*******************************************************************************
# Figure 3
#   Effects on Percent Voting by Mail and Turnout, Los Angeles County
#*******************************************************************************

# Left panel (Percent voting by mail)
Figure_PercVBM_by_UVBM_PAV_LA <- gen_figure(PercVBM_NPAV_Part1_est,
                                            PercVBM_PAV_Part1_est,
                                            y_axis_limits = c(0, 1),
                                            y_axis_breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1),
                                            y_axis_text = "% VBM",
                                            title_text = "")

# Right panel (Turnout)
Figure_Turnout_by_UVBM_PAV_LA <- gen_figure(Turnout_NPAV_Part1_est,
                                            Turnout_PAV_Part1_est,
                                            y_axis_limits = c(0.3, 0.5),
                                            y_axis_breaks = c(0.3, 0.35, 0.4, 0.45, 0.5),
                                            y_axis_text = "Turnout",
                                            title_text = "")

# Full Figure
Figure_combined_by_UVBM_PAV_LA <-
  plot_spacer() + Figure_PercVBM_by_UVBM_PAV_LA + plot_spacer() + Figure_Turnout_by_UVBM_PAV_LA + plot_spacer() +
  plot_layout(widths = c(1, 3, 1, 3, 1)) +
  plot_layout(guides = "collect") & theme(legend.position = "bottom", legend.direction = "vertical")

# Save Figure
ggsave("./Replication Package/Figure 3.pdf", Figure_combined_by_UVBM_PAV_LA,
       width = 12, height = 8, units = "in", dpi = 600)

#*******************************************************************************
# Table SI.5
#   Effects on Percent Voting by Mail, Los Angeles County
#*******************************************************************************

# Top Panel (Non-Permanent VBM Voters)
Table_SI5_NPAV <- data.frame(
  Voter_Type = c("Non-Permanent VBM Voters", NA_character_, NA_character_, NA_character_),
  N_Voter_Type = c(N_PercVBM_NPAV, NA_real_, NA_real_, NA_real_),
  Dist_Type = c("Universal VBM Districts", NA_character_, "Non-Universal VBM Districts", NA_character_),
  Level_2016 = c(PercVBM_NPAV_Part1_est$coef[2], NA_real_, PercVBM_NPAV_Part1_est$coef[1], NA_real_),
  Level_2020 = c(PercVBM_NPAV_Part1_est$coef[4], NA_real_, PercVBM_NPAV_Part1_est$coef[3], NA_real_),
  Difference = c(PercVBM_NPAV_Part2_est$coef[4], PercVBM_NPAV_Part2_est$se[4],
                 PercVBM_NPAV_Part2_est$coef[3], PercVBM_NPAV_Part2_est$se[3]),
  Diff_in_Diffs = c(PercVBM_NPAV_Part3_est$coef[4], PercVBM_NPAV_Part3_est$se[4], NA_real_, NA_real_)
)

# Middle Panel (Permanent VBM Voters)
Table_SI5_PAV <- data.frame(
  Voter_Type = c("Permanent VBM Voters", NA_character_, NA_character_, NA_character_),
  N_Voter_Type = c(N_PercVBM_PAV, NA_real_, NA_real_, NA_real_),
  Dist_Type = c("Universal VBM Districts", NA_character_, "Non-Universal VBM Districts", NA_character_),
  Level_2016 = c(PercVBM_PAV_Part1_est$coef[2], NA_real_, PercVBM_PAV_Part1_est$coef[1], NA_real_),
  Level_2020 = c(PercVBM_PAV_Part1_est$coef[4], NA_real_, PercVBM_PAV_Part1_est$coef[3], NA_real_),
  Difference = c(PercVBM_PAV_Part2_est$coef[4], PercVBM_PAV_Part2_est$se[4],
                 PercVBM_PAV_Part2_est$coef[3], PercVBM_PAV_Part2_est$se[3]),
  Diff_in_Diffs = c(PercVBM_PAV_Part3_est$coef[4], PercVBM_PAV_Part3_est$se[4], NA_real_, NA_real_)
)

# Bottom Panel (Non-Perm. And Perm. VBM Voters)
Table_SI5_ALL <- data.frame(
  Voter_Type = c("Non-Perm. And Perm. VBM Voters", NA_character_, NA_character_, NA_character_),
  N_Voter_Type = c(N_PercVBM_ALL, NA_real_, NA_real_, NA_real_),
  Dist_Type = c("Universal VBM Districts", NA_character_, "Non-Universal VBM Districts", NA_character_),
  Level_2016 = c(PercVBM_ALL_Part1_est$coef[2], NA_real_, PercVBM_ALL_Part1_est$coef[1], NA_real_),
  Level_2020 = c(PercVBM_ALL_Part1_est$coef[4], NA_real_, PercVBM_ALL_Part1_est$coef[3], NA_real_),
  Difference = c(PercVBM_ALL_Part2_est$coef[4], PercVBM_ALL_Part2_est$se[4],
                 PercVBM_ALL_Part2_est$coef[3], PercVBM_ALL_Part2_est$se[3]),
  Diff_in_Diffs = c(PercVBM_ALL_Part3_est$coef[4], PercVBM_ALL_Part3_est$se[4], NA_real_, NA_real_)
)

# Full Table
Table_SI5 <- bind_rows(Table_SI5_NPAV, Table_SI5_PAV, Table_SI5_ALL) %>%
  mutate(Level_2016 = round(100*Level_2016, 1),
         Level_2020 = round(100*Level_2020, 1),
         Difference = round(100*Difference, 2),
         Diff_in_Diffs = round(100*Diff_in_Diffs, 2))

# Save Table
write.csv(Table_SI5, "./Replication Package/Table SI.5.csv", row.names = FALSE, na = '')

#*******************************************************************************
# Table SI.6
#   Effects on Turnout, Los Angeles County
#*******************************************************************************

# Top Panel (Non-Permanent VBM Voters)
Table_SI6_NPAV <- data.frame(
  Voter_Type = c("Non-Permanent VBM Voters", NA_character_, NA_character_, NA_character_),
  N_Voter_Type = c(N_Turnout_NPAV, NA_real_, NA_real_, NA_real_),
  Dist_Type = c("Universal VBM Districts", NA_character_, "Non-Universal VBM Districts", NA_character_),
  Level_2016 = c(Turnout_NPAV_Part1_est$coef[2], NA_real_, Turnout_NPAV_Part1_est$coef[1], NA_real_),
  Level_2020 = c(Turnout_NPAV_Part1_est$coef[4], NA_real_, Turnout_NPAV_Part1_est$coef[3], NA_real_),
  Difference = c(Turnout_NPAV_Part2_est$coef[4], Turnout_NPAV_Part2_est$se[4],
                 Turnout_NPAV_Part2_est$coef[3], Turnout_NPAV_Part2_est$se[3]),
  Diff_in_Diffs = c(Turnout_NPAV_Part3_est$coef[4], Turnout_NPAV_Part3_est$se[4], NA_real_, NA_real_)
)

# Middle Panel (Permanent VBM Voters)
Table_SI6_PAV <- data.frame(
  Voter_Type = c("Permanent VBM Voters", NA_character_, NA_character_, NA_character_),
  N_Voter_Type = c(N_Turnout_PAV, NA_real_, NA_real_, NA_real_),
  Dist_Type = c("Universal VBM Districts", NA_character_, "Non-Universal VBM Districts", NA_character_),
  Level_2016 = c(Turnout_PAV_Part1_est$coef[2], NA_real_, Turnout_PAV_Part1_est$coef[1], NA_real_),
  Level_2020 = c(Turnout_PAV_Part1_est$coef[4], NA_real_, Turnout_PAV_Part1_est$coef[3], NA_real_),
  Difference = c(Turnout_PAV_Part2_est$coef[4], Turnout_PAV_Part2_est$se[4],
                 Turnout_PAV_Part2_est$coef[3], Turnout_PAV_Part2_est$se[3]),
  Diff_in_Diffs = c(Turnout_PAV_Part3_est$coef[4], Turnout_PAV_Part3_est$se[4], NA_real_, NA_real_)
)

# Bottom Panel (Non-Perm. And Perm. VBM Voters)
Table_SI6_ALL <- data.frame(
  Voter_Type = c("Non-Perm. And Perm. VBM Voters", NA_character_, NA_character_, NA_character_),
  N_Voter_Type = c(N_Turnout_ALL, NA_real_, NA_real_, NA_real_),
  Dist_Type = c("Universal VBM Districts", NA_character_, "Non-Universal VBM Districts", NA_character_),
  Level_2016 = c(Turnout_ALL_Part1_est$coef[2], NA_real_, Turnout_ALL_Part1_est$coef[1], NA_real_),
  Level_2020 = c(Turnout_ALL_Part1_est$coef[4], NA_real_, Turnout_ALL_Part1_est$coef[3], NA_real_),
  Difference = c(Turnout_ALL_Part2_est$coef[4], Turnout_ALL_Part2_est$se[4],
                 Turnout_ALL_Part2_est$coef[3], Turnout_ALL_Part2_est$se[3]),
  Diff_in_Diffs = c(Turnout_ALL_Part3_est$coef[4], Turnout_ALL_Part3_est$se[4], NA_real_, NA_real_)
)

# Full Table
Table_SI6 <- bind_rows(Table_SI6_NPAV, Table_SI6_PAV, Table_SI6_ALL) %>%
  mutate(Level_2016 = round(100*Level_2016, 1),
         Level_2020 = round(100*Level_2020, 1),
         Difference = round(100*Difference, 2),
         Diff_in_Diffs = round(100*Diff_in_Diffs, 2))

# Save Table
write.csv(Table_SI6, "./Replication Package/Table SI.6.csv", row.names = FALSE, na = '')

#*******************************************************************************
# Difference-in-Differences Analysis without Covariates by Party
# Results Reported in Figure 5, Table SI.8
#*******************************************************************************

# Registered Democrats
# Non-Permanent VBM Voters
Turnout_NPAV_DEM_Part1 <- run_model("Turnout", "Levels", "Non-Perm. VBM", "Democratic")
Turnout_NPAV_DEM_Part2 <- run_model("Turnout", "Differences", "Non-Perm. VBM", "Democratic")
Turnout_NPAV_DEM_Part3 <- run_model("Turnout", "Diff-in-Diffs", "Non-Perm. VBM", "Democratic")

# Permanent VBM Voters
Turnout_PAV_DEM_Part1 <- run_model("Turnout", "Levels", "Perm. VBM", "Democratic")
Turnout_PAV_DEM_Part2 <- run_model("Turnout", "Differences", "Perm. VBM", "Democratic")
Turnout_PAV_DEM_Part3 <- run_model("Turnout", "Diff-in-Diffs", "Perm. VBM", "Democratic")

# Non-Perm. and Perm. VBM Voters
Turnout_ALL_DEM_Part1 <- run_model("Turnout", "Levels", "All Voters", "Democratic")
Turnout_ALL_DEM_Part2 <- run_model("Turnout", "Differences", "All Voters", "Democratic")
Turnout_ALL_DEM_Part3 <- run_model("Turnout", "Diff-in-Diffs", "All Voters", "Democratic")

# Obtain Estimates and Clustered Standard Errors
# Non-Permanent VBM Voters
Turnout_NPAV_DEM_Part1_est <- get_estimates(Turnout_NPAV_DEM_Part1)
Turnout_NPAV_DEM_Part2_est <- get_estimates(Turnout_NPAV_DEM_Part2)
Turnout_NPAV_DEM_Part3_est <- get_estimates(Turnout_NPAV_DEM_Part3)

# Permanent VBM Voters
Turnout_PAV_DEM_Part1_est <- get_estimates(Turnout_PAV_DEM_Part1)
Turnout_PAV_DEM_Part2_est <- get_estimates(Turnout_PAV_DEM_Part2)
Turnout_PAV_DEM_Part3_est <- get_estimates(Turnout_PAV_DEM_Part3)

# Non-Perm. and Perm. VBM Voters
Turnout_ALL_DEM_Part1_est <- get_estimates(Turnout_ALL_DEM_Part1)
Turnout_ALL_DEM_Part2_est <- get_estimates(Turnout_ALL_DEM_Part2)
Turnout_ALL_DEM_Part3_est <- get_estimates(Turnout_ALL_DEM_Part3)

# Number of Observations
N_Turnout_NPAV_DEM <- nobs(Turnout_NPAV_DEM_Part1)
N_Turnout_PAV_DEM <- nobs(Turnout_PAV_DEM_Part1)
N_Turnout_ALL_DEM <- nobs(Turnout_ALL_DEM_Part1)

# Registered Republicans
# Non-Permanent VBM Voters
Turnout_NPAV_REP_Part1 <- run_model("Turnout", "Levels", "Non-Perm. VBM", "Republican")
Turnout_NPAV_REP_Part2 <- run_model("Turnout", "Differences", "Non-Perm. VBM", "Republican")
Turnout_NPAV_REP_Part3 <- run_model("Turnout", "Diff-in-Diffs", "Non-Perm. VBM", "Republican")

# Permanent VBM Voters
Turnout_PAV_REP_Part1 <- run_model("Turnout", "Levels", "Perm. VBM", "Republican")
Turnout_PAV_REP_Part2 <- run_model("Turnout", "Differences", "Perm. VBM", "Republican")
Turnout_PAV_REP_Part3 <- run_model("Turnout", "Diff-in-Diffs", "Perm. VBM", "Republican")

# Non-Perm. and Perm. VBM Voters
Turnout_ALL_REP_Part1 <- run_model("Turnout", "Levels", "All Voters", "Republican")
Turnout_ALL_REP_Part2 <- run_model("Turnout", "Differences", "All Voters", "Republican")
Turnout_ALL_REP_Part3 <- run_model("Turnout", "Diff-in-Diffs", "All Voters", "Republican")

# Obtain Estimates and Clustered Standard Errors
# Non-Permanent VBM Voters
Turnout_NPAV_REP_Part1_est <- get_estimates(Turnout_NPAV_REP_Part1)
Turnout_NPAV_REP_Part2_est <- get_estimates(Turnout_NPAV_REP_Part2)
Turnout_NPAV_REP_Part3_est <- get_estimates(Turnout_NPAV_REP_Part3)

# Permanent VBM Voters
Turnout_PAV_REP_Part1_est <- get_estimates(Turnout_PAV_REP_Part1)
Turnout_PAV_REP_Part2_est <- get_estimates(Turnout_PAV_REP_Part2)
Turnout_PAV_REP_Part3_est <- get_estimates(Turnout_PAV_REP_Part3)

# Non-Perm. and Perm. VBM Voters
Turnout_ALL_REP_Part1_est <- get_estimates(Turnout_ALL_REP_Part1)
Turnout_ALL_REP_Part2_est <- get_estimates(Turnout_ALL_REP_Part2)
Turnout_ALL_REP_Part3_est <- get_estimates(Turnout_ALL_REP_Part3)

# Number of Observations
N_Turnout_NPAV_REP <- nobs(Turnout_NPAV_REP_Part1)
N_Turnout_PAV_REP <- nobs(Turnout_PAV_REP_Part1)
N_Turnout_ALL_REP <- nobs(Turnout_ALL_REP_Part1)

# Registered NPA Voters
# Non-Permanent VBM Voters
Turnout_NPAV_NPA_Part1 <- run_model("Turnout", "Levels", "Non-Perm. VBM", "NPA")
Turnout_NPAV_NPA_Part2 <- run_model("Turnout", "Differences", "Non-Perm. VBM", "NPA")
Turnout_NPAV_NPA_Part3 <- run_model("Turnout", "Diff-in-Diffs", "Non-Perm. VBM", "NPA")

# Permanent VBM Voters
Turnout_PAV_NPA_Part1 <- run_model("Turnout", "Levels", "Perm. VBM", "NPA")
Turnout_PAV_NPA_Part2 <- run_model("Turnout", "Differences", "Perm. VBM", "NPA")
Turnout_PAV_NPA_Part3 <- run_model("Turnout", "Diff-in-Diffs", "Perm. VBM", "NPA")

# Non-Perm. and Perm. VBM Voters
Turnout_ALL_NPA_Part1 <- run_model("Turnout", "Levels", "All Voters", "NPA")
Turnout_ALL_NPA_Part2 <- run_model("Turnout", "Differences", "All Voters", "NPA")
Turnout_ALL_NPA_Part3 <- run_model("Turnout", "Diff-in-Diffs", "All Voters", "NPA")

# Obtain Estimates and Clustered Standard Errors
# Non-Permanent VBM Voters
Turnout_NPAV_NPA_Part1_est <- get_estimates(Turnout_NPAV_NPA_Part1)
Turnout_NPAV_NPA_Part2_est <- get_estimates(Turnout_NPAV_NPA_Part2)
Turnout_NPAV_NPA_Part3_est <- get_estimates(Turnout_NPAV_NPA_Part3)

# Permanent VBM Voters
Turnout_PAV_NPA_Part1_est <- get_estimates(Turnout_PAV_NPA_Part1)
Turnout_PAV_NPA_Part2_est <- get_estimates(Turnout_PAV_NPA_Part2)
Turnout_PAV_NPA_Part3_est <- get_estimates(Turnout_PAV_NPA_Part3)

# Non-Perm. and Perm. VBM Voters
Turnout_ALL_NPA_Part1_est <- get_estimates(Turnout_ALL_NPA_Part1)
Turnout_ALL_NPA_Part2_est <- get_estimates(Turnout_ALL_NPA_Part2)
Turnout_ALL_NPA_Part3_est <- get_estimates(Turnout_ALL_NPA_Part3)

# Number of Observations
N_Turnout_NPAV_NPA <- nobs(Turnout_NPAV_NPA_Part1)
N_Turnout_PAV_NPA <- nobs(Turnout_PAV_NPA_Part1)
N_Turnout_ALL_NPA <- nobs(Turnout_ALL_NPA_Part1)

#*******************************************************************************
# Figure 5
#   Effects on Turnout by Party Registration, Los Angeles County
#*******************************************************************************

# Left Panel (Registered Democrats)
Figure_Turnout_DEM_by_UVBM_PAV_LA <- gen_figure(Turnout_NPAV_DEM_Part1_est,
                                                Turnout_PAV_DEM_Part1_est,
                                                y_axis_limits = c(0.2, 0.56),
                                                y_axis_breaks = seq(0.2, 0.55, by = 0.05),
                                                y_axis_text = "Turnout",
                                                title_text = "DEM Voters")

# Middle Panel (Registered Republicans)
Figure_Turnout_REP_by_UVBM_PAV_LA <- gen_figure(Turnout_NPAV_REP_Part1_est,
                                                Turnout_PAV_REP_Part1_est,
                                                y_axis_limits = c(0.2, 0.56),
                                                y_axis_breaks = seq(0.2, 0.55, by = 0.05),
                                                y_axis_text = "Turnout",
                                                title_text = "REP Voters")

# Right Panel (Registered NPA Voters)
Figure_Turnout_NPA_by_UVBM_PAV_LA <- gen_figure(Turnout_NPAV_NPA_Part1_est,
                                                Turnout_PAV_NPA_Part1_est,
                                                y_axis_limits = c(0.2, 0.56),
                                                y_axis_breaks = seq(0.2, 0.55, by = 0.05),
                                                y_axis_text = "Turnout",
                                                title_text = "NPA Voters")

# Full Figure
Figure_Turnout_Parties_by_UVBM_PAV_LA <-
  Figure_Turnout_DEM_by_UVBM_PAV_LA + Figure_Turnout_REP_by_UVBM_PAV_LA + Figure_Turnout_NPA_by_UVBM_PAV_LA +
  plot_layout(guides = "collect") & theme(legend.position = "bottom", legend.direction = "vertical")

# Save Figure
ggsave("./Replication Package/Figure 5.pdf",
       Figure_Turnout_Parties_by_UVBM_PAV_LA, width = 12, height = 8, units = "in", dpi = 600)

#*******************************************************************************
# Table SI.8
#   Effects on Turnout by Party Registration, Los Angeles County
#*******************************************************************************

# Top Left Panel (Registered Democrats, Non-Permanent VBM Voters)
Table_SI8_NPAV_DEM <- data.frame(
  Voter_Type = c("Non-Permanent VBM Voters", NA_character_, NA_character_, NA_character_),
  N_Voter_Type = c(N_Turnout_NPAV_DEM, NA_real_, NA_real_, NA_real_),
  Dist_Type = c("Universal VBM Districts", NA_character_, "Non-Universal VBM Districts", NA_character_),
  Level_2016 = c(Turnout_NPAV_DEM_Part1_est$coef[2], NA_real_, Turnout_NPAV_DEM_Part1_est$coef[1], NA_real_),
  Level_2020 = c(Turnout_NPAV_DEM_Part1_est$coef[4], NA_real_, Turnout_NPAV_DEM_Part1_est$coef[3], NA_real_),
  Diff_in_Diffs = c(Turnout_NPAV_DEM_Part3_est$coef[4], Turnout_NPAV_DEM_Part3_est$se[4], NA_real_, NA_real_)
)

# Center Left Panel (Registered Democrats, Permanent VBM Voters)
Table_SI8_PAV_DEM <- data.frame(
  Voter_Type = c("Permanent VBM Voters", NA_character_, NA_character_, NA_character_),
  N_Voter_Type = c(N_Turnout_PAV_DEM, NA_real_, NA_real_, NA_real_),
  Dist_Type = c("Universal VBM Districts", NA_character_, "Non-Universal VBM Districts", NA_character_),
  Level_2016 = c(Turnout_PAV_DEM_Part1_est$coef[2], NA_real_, Turnout_PAV_DEM_Part1_est$coef[1], NA_real_),
  Level_2020 = c(Turnout_PAV_DEM_Part1_est$coef[4], NA_real_, Turnout_PAV_DEM_Part1_est$coef[3], NA_real_),
  Diff_in_Diffs = c(Turnout_PAV_DEM_Part3_est$coef[4], Turnout_PAV_DEM_Part3_est$se[4], NA_real_, NA_real_)
)

# Bottom Left Panel (Registered Democrats, Non-Perm. and Perm. VBM Voters)
Table_SI8_ALL_DEM <- data.frame(
  Voter_Type = c("Non-Perm. And Perm. VBM Voters", NA_character_, NA_character_, NA_character_),
  N_Voter_Type = c(N_Turnout_ALL_DEM, NA_real_, NA_real_, NA_real_),
  Dist_Type = c("Universal VBM Districts", NA_character_, "Non-Universal VBM Districts", NA_character_),
  Level_2016 = c(Turnout_ALL_DEM_Part1_est$coef[2], NA_real_, Turnout_ALL_DEM_Part1_est$coef[1], NA_real_),
  Level_2020 = c(Turnout_ALL_DEM_Part1_est$coef[4], NA_real_, Turnout_ALL_DEM_Part1_est$coef[3], NA_real_),
  Diff_in_Diffs = c(Turnout_ALL_DEM_Part3_est$coef[4], Turnout_ALL_DEM_Part3_est$se[4], NA_real_, NA_real_)
)

# Full Left Panel (Registered Democrats)
Table_SI8_DEM <- bind_rows(Table_SI8_NPAV_DEM, Table_SI8_PAV_DEM, Table_SI8_ALL_DEM) %>%
  mutate(Level_2016 = round(100*Level_2016, 1),
         Level_2020 = round(100*Level_2020, 1),
         Diff_in_Diffs = round(100*Diff_in_Diffs, 2))

# Top Middle Panel (Registered Republicans, Non-Permanent VBM Voters)
Table_SI8_NPAV_REP <- data.frame(
  Voter_Type = c("Non-Permanent VBM Voters", NA_character_, NA_character_, NA_character_),
  N_Voter_Type = c(N_Turnout_NPAV_REP, NA_real_, NA_real_, NA_real_),
  Dist_Type = c("Universal VBM Districts", NA_character_, "Non-Universal VBM Districts", NA_character_),
  Level_2016 = c(Turnout_NPAV_REP_Part1_est$coef[2], NA_real_, Turnout_NPAV_REP_Part1_est$coef[1], NA_real_),
  Level_2020 = c(Turnout_NPAV_REP_Part1_est$coef[4], NA_real_, Turnout_NPAV_REP_Part1_est$coef[3], NA_real_),
  Diff_in_Diffs = c(Turnout_NPAV_REP_Part3_est$coef[4], Turnout_NPAV_REP_Part3_est$se[4], NA_real_, NA_real_)
)

# Center Middle Panel (Registered Republicans, Permanent VBM Voters)
Table_SI8_PAV_REP <- data.frame(
  Voter_Type = c("Permanent VBM Voters", NA_character_, NA_character_, NA_character_),
  N_Voter_Type = c(N_Turnout_PAV_REP, NA_real_, NA_real_, NA_real_),
  Dist_Type = c("Universal VBM Districts", NA_character_, "Non-Universal VBM Districts", NA_character_),
  Level_2016 = c(Turnout_PAV_REP_Part1_est$coef[2], NA_real_, Turnout_PAV_REP_Part1_est$coef[1], NA_real_),
  Level_2020 = c(Turnout_PAV_REP_Part1_est$coef[4], NA_real_, Turnout_PAV_REP_Part1_est$coef[3], NA_real_),
  Diff_in_Diffs = c(Turnout_PAV_REP_Part3_est$coef[4], Turnout_PAV_REP_Part3_est$se[4], NA_real_, NA_real_)
)

# Bottom Middle Panel (Registered Republicans, Non-Perm. and Perm. VBM Voters)
Table_SI8_ALL_REP <- data.frame(
  Voter_Type = c("Non-Perm. And Perm. VBM Voters", NA_character_, NA_character_, NA_character_),
  N_Voter_Type = c(N_Turnout_ALL_REP, NA_real_, NA_real_, NA_real_),
  Dist_Type = c("Universal VBM Districts", NA_character_, "Non-Universal VBM Districts", NA_character_),
  Level_2016 = c(Turnout_ALL_REP_Part1_est$coef[2], NA_real_, Turnout_ALL_REP_Part1_est$coef[1], NA_real_),
  Level_2020 = c(Turnout_ALL_REP_Part1_est$coef[4], NA_real_, Turnout_ALL_REP_Part1_est$coef[3], NA_real_),
  Diff_in_Diffs = c(Turnout_ALL_REP_Part3_est$coef[4], Turnout_ALL_REP_Part3_est$se[4], NA_real_, NA_real_)
)

# Full Middle Panel (Registered Republicans)
Table_SI8_REP <- bind_rows(Table_SI8_NPAV_REP, Table_SI8_PAV_REP, Table_SI8_ALL_REP) %>%
  mutate(Level_2016 = round(100*Level_2016, 1),
         Level_2020 = round(100*Level_2020, 1),
         Diff_in_Diffs = round(100*Diff_in_Diffs, 2))

# Top Right Panel (Registered NPA Voters, Non-Permanent VBM Voters)
Table_SI8_NPAV_NPA <- data.frame(
  Voter_Type = c("Non-Permanent VBM Voters", NA_character_, NA_character_, NA_character_),
  N_Voter_Type = c(N_Turnout_NPAV_NPA, NA_real_, NA_real_, NA_real_),
  Dist_Type = c("Universal VBM Districts", NA_character_, "Non-Universal VBM Districts", NA_character_),
  Level_2016 = c(Turnout_NPAV_NPA_Part1_est$coef[2], NA_real_, Turnout_NPAV_NPA_Part1_est$coef[1], NA_real_),
  Level_2020 = c(Turnout_NPAV_NPA_Part1_est$coef[4], NA_real_, Turnout_NPAV_NPA_Part1_est$coef[3], NA_real_),
  Diff_in_Diffs = c(Turnout_NPAV_NPA_Part3_est$coef[4], Turnout_NPAV_NPA_Part3_est$se[4], NA_real_, NA_real_)
)

# Center Right Panel (Registered NPA Voters, Permanent VBM Voters)
Table_SI8_PAV_NPA <- data.frame(
  Voter_Type = c("Permanent VBM Voters", NA_character_, NA_character_, NA_character_),
  N_Voter_Type = c(N_Turnout_PAV_NPA, NA_real_, NA_real_, NA_real_),
  Dist_Type = c("Universal VBM Districts", NA_character_, "Non-Universal VBM Districts", NA_character_),
  Level_2016 = c(Turnout_PAV_NPA_Part1_est$coef[2], NA_real_, Turnout_PAV_NPA_Part1_est$coef[1], NA_real_),
  Level_2020 = c(Turnout_PAV_NPA_Part1_est$coef[4], NA_real_, Turnout_PAV_NPA_Part1_est$coef[3], NA_real_),
  Diff_in_Diffs = c(Turnout_PAV_NPA_Part3_est$coef[4], Turnout_PAV_NPA_Part3_est$se[4], NA_real_, NA_real_)
)

# Bottom Right Panel (Registered NPA Voters, Non-Perm. and Perm. VBM Voters)
Table_SI8_ALL_NPA <- data.frame(
  Voter_Type = c("Non-Perm. And Perm. VBM Voters", NA_character_, NA_character_, NA_character_),
  N_Voter_Type = c(N_Turnout_ALL_NPA, NA_real_, NA_real_, NA_real_),
  Dist_Type = c("Universal VBM Districts", NA_character_, "Non-Universal VBM Districts", NA_character_),
  Level_2016 = c(Turnout_ALL_NPA_Part1_est$coef[2], NA_real_, Turnout_ALL_NPA_Part1_est$coef[1], NA_real_),
  Level_2020 = c(Turnout_ALL_NPA_Part1_est$coef[4], NA_real_, Turnout_ALL_NPA_Part1_est$coef[3], NA_real_),
  Diff_in_Diffs = c(Turnout_ALL_NPA_Part3_est$coef[4], Turnout_ALL_NPA_Part3_est$se[4], NA_real_, NA_real_)
)

# Full Right Panel (Registered NPA Voters)
Table_SI8_NPA <- bind_rows(Table_SI8_NPAV_NPA, Table_SI8_PAV_NPA, Table_SI8_ALL_NPA) %>%
  mutate(Level_2016 = round(100*Level_2016, 1),
         Level_2020 = round(100*Level_2020, 1),
         Diff_in_Diffs = round(100*Diff_in_Diffs, 2))

# Full Table
Table_SI8 <- bind_cols(Table_SI8_DEM, Table_SI8_REP, Table_SI8_NPA)

# Save Table
write.csv(Table_SI8, file = "./Replication Package/Table SI.8.csv", row.names = FALSE, na = '')
